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Abstract 

Highly accurate nonrelativistic ground-state wave function and energy of the lithium atom is obtained 
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evaluation of Hylleraas integrals with the help of recursion relations. 
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I. INTRODUCTION 



Theoretical predictions for the energy levels of light few-electron atoms are much less accurate 
than for the hydrogenic systems. It is for two reasons. The nonrelativistic wave function has to in- 
clude electron correlations to a high degree of accuracy. This can be achieved by using a Hylleraas 
basis set, but it is quite difficult to evaluate integrals with Hylleraas functions for three and more 
electrons. The second reason is the difficulty in the accurate treatment of relativistic and radiative 
corrections. The commonly used Dirac-Coulomb Hamiltonian for few-electron atoms does not 
include relativistic corrections properly as it cannot be derived from quantum electrodynamic the- 
ory and its continuous spectrum ranges from — oo to +00. One of the possible approaches is the 
derivation of an effective Hamiltonian [ 1 ] within the so called NRQED theory. Matrix elements of 
this Hamiltonian give exact correction to the energy at specified order in the fine structure constant 
a. However, this Hamiltonian becomes quite complicated at higher orders and for example m a 6 
corrections has been obtained for few low lying states of helium only J21 0], not for lithium nor 
beryllium atoms. 

Theoretical predictions for light hydrogen-like atoms are at present limited by uncertainty in 
higher-order two-loop electron self-energy corrections [4], which is a few kHz for the IS state. For 
helium-like atoms predictions are approximately 10 3 times less accurate. Since, the nonrelativistic 
wave function was computed very accurately using Hylleraas [Q] or exponential basis sets Q], 
the uncertainty in its levels comes mainly from the unknown ma 7 terms. These corrections are 
currently under investigation in the context of helium 2 3 Pj fine splitting. For lithium atoms, 
the Hylleraas functions give very accurate nonrelativistic wave function and energies [7], but the 
precise calculation of three-electron integrals with Hylleraas functions is very time consuming 
IE |2D, and so far no result for ma 6 corrections have been obtained. For the beryllium atom 



the most accurate results have been obtained with explicitly correlated Gaussian functions [10]. 
Although it was possible to calculate accurately the leading relativistic and QED corrections 111 ill , 
the final accuracy is limited by the nonrelativistic energy. Moreover, this basis cannot be used for 
higher order corrections since Gaussian wave functions do not fulfill the cusp condition. 

So far the most accurate results for various states of the lithium atom were obtained by Yan and 
Drake in Ref. |Q]. Here, we present even more accurate results for the lithium ground state, as a 
demonstration of an analytic method to compute the integrals with Hylleraas functions 

O- This 

new method is based on recursion relations between integrals with different powers of electron- 
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nucleus and inter-electron distances, which are fast and numerically stable for generating large 
basis sets. Our result for the ground state energy 

E = -7.478 060 323 904 l(t$) , (1) 

is significantly below the previous one, obtained in which is E = —7.478 060 323 650 3(71). 
As a further application of the analytic approach, we obtain the leading relativistic corrections to 
the binding energy by the calculation of the expectation value of Breit-Pauli Hamiltonian in Eq. 
(fT3l . For this we used recursion relations for extended Hylleraas integrals with 1/r?- and 1/rf 
terms. They have been derived in Jl^l and in this work respectively. 

n 

In the next Section we construct the nonrelativistic wave function, similarly to Ref. [7] and 
obtain the ground state nonrelativistic energy and the wave function. In Sec. Ill we compute 
the leading relativistic correction as given by the Breit-Pauli Hamiltonian. In Sec. IV we derive 
recursion relations for Hylleraas integrals containing 1/rf which among others, are necessary 
for relativistic matrix elements. In Sec. V we summarize our result and present prospects for 
calculation of higher order terms as well as the calculation of Hylleraas integrals for 4 and more 
electrons. 



II. NONRELATIVISTIC WAVE FUNCTION AND ENERGY 

In the construction of the wave function we closely follow the works of Yan and Drake in 
The ground state wave function ^ is expressed as a linear combination of ip, the anti symmetrized 
product of and the spin function x 

ip = ^4[0(fi,f 2 ,f 3 ) x) , (2) 
0(ri, f 2 , f 3 ) = e~ wi ri ' W2 " 3 ra r n z l r™| r™ 4 r™ 5 r™ 6 , (3) 
X = a(l) (3(2) a(3) - (3(1) a(2) a(3) , (4) 

with all rij nonnegative integers and Wi E R + . The matrix element of the Hamiltonian H 

a=l x ' a>b=l 

or of any spin independent operator can be expressed after eliminating spin variables, as 

(mW) = (2 0(1, 2, 3) + 2 0(2, 1, 3) - 0(3, 1, 2) - 0(2, 3, 1) - 0(1, 3, 2) - 0(3, 2, 1) | 

# |0'(1,2,3)>. (6) 



In this way the calculation of this matrix elements is brought to Hylleraas integrals, namely the 
integrals with respect to of the form 

oPri f d 3 r 2 f d 3 r 3 



/(ni,n 2 ,n 3 ,n4,n 5 ,n 6 ) 



4:71 / 47T 



-wi r\— W2 T2— U13 rg 



Air 



~ni—l „n<z—\ „ns-\ „U4-1 „ns-l „ne-l 
'23 '31 '12 '1 '2 '3 



with nonnegative integers rii. These are performed analytical 



recursion relations for larger using formulas derived in II 1211 



ly for m, ri2, n 3 = 



1 



(7) 



and by 



The total wave function is generated from all <fi in Eq. © with rii satisfying condition 



^2 Hi - n ' 

i=l 



(8) 



for f2 between 3 and 12. For each we minimize energy with respect to the free parameters Wi in 
Eq. ©. We noticed that the use of only one set of u>j's does not lead to accurate results, therefore, 
following Yan and Drake [7], we divide the whole basis set into 5 sectors, each one with its own 
set of Wi's. This division goes as follows [7] 

sector 1: all n 3 , n\ = 0, n 2 = 0; 
sector 2: all n 3 , n\ = 0, n 2 7^ 0; 
sector 3: all n 3 , ri\ 7^ 0, n 2 = 0; 
sector 4: n 3 = 0, n\ 7^ 0, n 2 7^ 0; 
sector 5: n 3 7^ 0, ni 7^ 0, n 2 7^ 0; 

To avoid numerical instabilities, within each sector we drop the terms with n 4 > n 5 (or n 4 < n 5 ) 
and for n 4 = n 5 drop terms with n\ > n 2 (or n\ < n 2 ). This division allows for a significant 
improvements of nonrelativistic energies by optimization of all five sets of w^s. These nonlinear 
parameters are obtained by Newton method of searching zeros using analytic derivatives 



aw 



H 



aw 



dw 



(9) 



In the numerical calculations, we use sextuple precision for recursion relations and quadruple 
precision for all other arithmetics to obtain the wave function and the energy up to Q = 12. The 
results obtained for ground state energies are presented in Table I. The penultimate row is a result 
of extrapolation to infinite length of the basis set, and the last raw are previous results of Yan and 
Drake The result for the nonrelativistic energy is significantly below the previous estimate |Q] 
and indicates that extrapolation to infinite basis length does not always give the right result. In the 



TABLE I: Ground state nonrelativistic energies and expectation values of Dirac ^-functions obtained using 
Drachman formulae for various basis length. 



Q No. of terms E(Q) E„^a) £«>6<* 3 (»a&) 



3 


50 


-7.477 981 524089 7 


13 


843 446 803 98 


t~\ r" A A ASA ■ 1 C\i^ 

0.544 164 351 92 


4 


120 


-7.478 052 334 642 2 


13 


842 288 641 67 


0.544 331 564 16 


5 


256 


-7.478 059 463 915 8 


13 


O A r\ C f\C\ A 'I A / O 

842 509 17463 


r\ c A A T"4 ^ OTA AC 

0.544 327 87045 


6 


512 


-7.478 060 208 663 7 


13 


842 637 966 67 


0.544 325 26063 


( 


918 


-7.478 060 310 362 9 


1 r\ 

13 


842 606 662 38 


f\ C A A T1 /I TOO OC 

0.544 324788 85 


8 


1589 


-7 478 060 320 507 6 


13 


842 608 24076 


544 324 697 02 


9 


2625 


-7.478 060 323 4501 


13 


842610098 57 


0.544 32462945 


10 


4172 


-7.478 060 323 775 


13 


842610698 67 


0.544 324 627 57 


11 


6412 


-7.478 060 323 8610 


13 


842 61077919 


0.544 324 63150 


12 


9576 


-7.478 060 323 889 7 


13 


842 61078106 


0.544 324 63205 


oo 


oo 


-7.478 060 323 904 1(±^) 


13 


842 610783 46(100) 


0.544 324 633 96(50) 


Refs. [7, 15] 


oo 


-7.478 060 323 650 3(71) 


13 


842 609 642(55) 


0.544 329 79(31) 



same Table we present results for the Dirac 5 functions, which also differs from previous results in 
I15I1 . We observe, the the number of significant digits for Dirac 5 is increased by using Drachman 
formulae Jlfjj . namely 



47r(V\5 3 {r ab )\y) = 2 (V 
4tt<^ \5 3 (r a )\ tf) = 4/ 



—(Ev-V) 

Tab 



1 



(£* - V) 



*) -2 Wv c ^ 



1 



(10) 
(11) 



where V is a total potential energy in Eq. ©. 
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III. LEADING RELATIVISTIC CORRECTION TO BINDING ENERGY 



The leading relativistic corrections to energy levels are given by the expectation values of the 
Breit-Pauli Hamiltonian 



Pa , TT Z a .. 



+ 



2 7r a 



in" 2m 2 
7ra 



5%r a ) + 



Am z r„ 



2 m 



r 



«6 



3 m 2 



a <y % n o{ I „„■„■ _ r 



4 m 2 r\ h 



'i " li I ^ij g ' ab \ ab 

r ab 



+ 



a 



4 m 2 r^ b 



2 (<7 ■ f ab xp b -a b - r ab x p a ) + (a b • r a6 x p b - a a ■ r ab x p a ) 



(12) 



For states with vanishing angular momentum L and spin S = 1/2, the expectation value is simpli- 
fied to the form 



£(4) = (y\ H W\y) = /J2 



->,\ 

Pa 



7T Z a 
2 m 2 



5 3 (r a ) 



2 m' 



Pa 



+ 

Tab 



i 1 
' ab' ab 



r 



ab 



(13) 



has already been obtained in works 1115 . 



17fl . Calculations of these matrix elements involves 



the usual Hylleraas integrals with all rij nonnegative and extended integrals, namely with one pa- 
rameter rij equal to —1. The direct numerical method to calculate these integrals was presented 

nn 

in BSII2D- Here we apply the analytic approach. Recursion relations for the case of ri\ or or 
equal to —1 have been obtained in II 1 311 . Hylleraas integrals involving n A or n 5 or n e equal to —1 
can in principle be obtained by the integration of the usual Hylleraas integral with respect to the 
corresponding parameter w-i [13]. However, some recursion relations may become unstable, for 
example in the case of n A = — 1 the recursion in n\ is numerically unstable for large w\. To avoid 
this problem we derive in the next section stable recursion relations for extended Hylleraas inte- 
grals with rii = — 1 for i = 4, 5, 6. Numerical results for matrix elements of the Breit Hamiltonian 
using these recursion relations, has been presented in Table I and II. One observes that the lowest 
convergence is for the —p 4 /8 term, and in spite of the differences for separate matrix elements, the 
total relativistic correction is in good agreement with the former result in 111 511 . 
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TABLE II: Matrix elements of the Breit-Pauli Hamiltonian in atomic units. 



\l 


v-^ 1 vv4 

Z^a ~~ 8 V a 


V- 1 v-n /S ij , r lb r3 nh\\-73 
T,a>b2^a{r ab +^fM 

nb 


0-C41 


3 


-78.587 286 69090 


-0.438 632 545 84 


-12.080670 336 80 


4 


-78.557 331 85961 


-0.436 096 58640 


-12.053 111 944 61 


•5 


-78.556 355 905 97 


-0.435 697 34491 


-12.050709 116 55 


6 


-78.556 71450343 


-0.435 616426 50 


-12.050 388 076 38 


7 


-78.556 195 780 85 


-0.435 602 362 02 


-12.050004 29451 


8 


-78.556 162 642 13 


-0.435 599 523 90 


-12.049 961 162 16 


9 


-78.556 137477 61 


-0.435 598 217 44 


-12.049 926 149 76 


10 


-78.556135 73401 


-0.435 598 047 58 


-12.049 921414 27 


11 


-78.556 131 596 34 


-0.435 597 963 57 


-12.049 916 800 81 


12 


-78.556128 63210 


-0.435 597 91050 


-12.049 913 77296 


oo 


-78.556112 88(200) 


-0.435 597 765(50) 


-12.049 897 86(200) 


Ref. 1151 


-78.556135 55(148) 


-0.435 598 001(137) 


-12.049 909 94(180) 



IV. RECURSION RELATIONS FOR THREE-ELECTRON EXTENDED HYLLERAAS INTE- 
GRAL WITH \jr\ 

In the former section we calculated relativistic corrections. For this we needed various extended 
Hylleraas integrals, among them, integrals with 1/rf, which are being derived here. To obtain 
recursion relations for three-electron Hylleraas integral in Eq. Q, one first considers the integral 

G 



G(mi,m 2 ,m3;m4,m 5 ,m 6 ) = — J d 3 h J d 3 k 2 J d 3 k 3 (k 2 + ul) mi {k\ + u\) m2 

(k 2 3 + ul)~ m3 (*4 + wlY^ + wlY^ (k 2 21 + wl)~ m \ (14) 

which is related to / by: /(0, 0, 0, 0, 0, 0) = G(l, 1, 1, 1, 1, l)\ Ul =u 2 =u 3 =o- The following 9 inte- 
gration by part identities are valid because the integral of the derivative of a function vanishing at 
infinity vanishes, 



.2\-l 



= id(i, j) = j d 3 h I d 3 k 2 j d 3 k 3 kj (k\ + u\ 
J J J d ki 

{k\ + uir 1 (ki + ui)-\ki 2 + wir 1 (k 2 13 + wjy 1 (& + wir 1 



(15) 
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where i,j = 1,2,3. The reduction of the scalar products from the numerator leads to the identities 
for the linear combination of the G functions. If any of the arguments is equal to 0, then G becomes 
a known two-electron Hylleraas type integral. These identities are used to derive various recursion 
relations. Here, we derive a set of recursions for the case when n 4 , n 5 or n 6 is equal to —1. Let 
us assume that n 4 = —1. The analytic expression for /(0, 0,0, — l,n 5 ,n 6 ) involves powers of 
w 2 — w 3 in the denominator which is not very convenient in high precision numerical calculations. 
Instead, we use recursions for /(0, 0, 0, 0, n 5 , n 6 ) and numerically integrate with respect to w±, 
namely 

POO 

/(0, 0,0,-1, n 5 , n 6 ) = / d Wl f (0,0,0, 0,n 5 ,n 6 ). (16) 

These recursions are derived as follows. We take id(i, i) with % — 1, 2, 3 and put u { = 0. Result- 
ing three equations are solved against three unknowns: (7(1,1,1,2,1,1), (7(1,1,1,1,2,1), and 
(7(1, 1,1, 1,1, 2). The solution for the last two G functions is the following 

(7(1, 1,1,1,2,1) = ^ [ G (°> !> X > !> !> 2 ) " ^ °> X > 2 ) " ^ °> 2 > X ) 

+<7(1, 1, 0, 2, 1, 1) + G(l, 1, 1, 1, 1, l)/2] , (17) 

G(l, 1,1,1,1,2) = ^ [G(0, 1, 1, 1, 2, 1) + G(l, 0, 1, 2, 1, 1) - G(l, 1, 0, 1, 2, 1) 

-(7(1,1,0, 2, l,l) + G(l,l,l,l,l,l)/2]. (18) 

By differentiation with respect to w 2 and w 3 one obtains the following recursion relations 

/(0, 0,0,0, n 5 + l,n 6 ) = [(n 5 + l)/(0, 0,0, 0,n 5 , ne)^!^ 

-(n 5 + 1) n 6 /(0, 0, 0, 0, n 5 , n 6 - 1) wi 

+n 6 /(0, 0, 0, 0, n 5 + 1, n 6 - 1) wi w 2 

-n 6 F(n 5 , n 6 -l, -1, u>i + w 2 , w 3 , 0) 

+n 6 T(n 6 - 1, 715, -1, wi + w 3 , w 2 , 0) 

-r(n 6 , n 5 , -1, wi + w 3 , w 2 , 0) wi 

+r(n 5 + 7i 6 , 0, -1, u> 2 + w 3 , 0) iui 

+r(n 5 , n 6 , -1, iui + tu 2 , w 3 , 0) w 3 

-r(n 6 ,7i 5 , -l,wi + w 3 ,w 2 ,0) w 3 ] , (19) 
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/(O, 0,0,0,715,7*6 + 1) = [(n 6 + l) f (0,0,0, 0,n 5 ,n 6 ) Wl w 2 

Wi w 2 w 3 L 

-n 5 (n 6 + 1) f(0, 0, 0, 0, n 5 - 1, n 6 ) w x 

+n 5 f(0, 0, 0, 0, n 5 - 1, n 6 + 1) w x w 3 

+n 5 T(n 5 - 1, 7i 6 , -1, Wi + w 2 , w 3 , 0) 

-rz 5 r(n 6 , n 5 - 1, -1, + ty 3 , ty 2 , 0) 

-r(n 5 , n 6 , -1, wi + U7 2 , u> 3 , 0) u>i 

+r(n 5 + n 6 , 0, -1, w 2 + w 3 , wi, 0) wi 

-r(n 5 , n 6 , -1, wi + w 2 , w 3 , 0) w 2 

+r(n 6 ,7i 5 , -l,wt + w 3 ,w 2 ,0) w 2 ] . (20) 

where T is a known two-electron integral 



r(n 1 ,n 2 ,7i3,a a ,a 2 ,a 3 ) = / ^ / ^^n-a.n-cam^-i^-i^-^ (21) 



The integration in Eq. dT6t is performed numerically using adapted points and weights to the 
function which has logarithmic end-point singularity, namely 

/ dx [Wi(x) + W 2 (x) ln(x)] , (22) 
Jo 

where Wi are functions without any singularities. The method to obtain n adapted points and 
weights is presented in Appendix A, and this integral is exact for Wi being polynomials up to the 
order n — 1. In the actual calculations we achieved at least 48 digits precision with only 100 points. 
Having obtained f(0, 0, 0, —1, n 5 , n 6 ) we construct recursion relations in n\, n 2 , and n 3 . This is 
achieved in two steps. In the first step we use integration by parts in momentum representation Eq. 
(fT5t . to form the following linear combination 

id(2, 2) + id(3, 3) - id(l, 1) = 2 [G(0, 1, 1, 1, 1, 2) + G(0, 1, 1, 1, 2, 1) - G(l, 0,1,1, 1, 2) 

-G(l, 1, 0, 1, 2, 1) - G(l, 1, 1, 1, 1, l)/2 - G(2, 1, 1, 1, 1, 1) u\ 
-G(l, 1, 1, 1, 1, 2) (uj - u\) + Gil, 2, 1, 1, 1, 1) u\ 
-Gil, 1, 1, 1, 2, 1) iu\ - ul) + G(l, 1, 2, 1, 1, 1) uj 
+£(1,1, 1,2,1,1)™?] =0. (23) 
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We integrate with respect to wi and differentiate over u±, u 2 , u 3 , w 2 , and w 3 to obtain the main 
formula 

f(n 1 ,n 2 ,n 3 , -l,n 5 ,n 6 ) = ■ — r [ 

\n 2 + n 3 -n 1 )w 2 w 3 

(n-y - 1) ni n 5 f(ni - 2, n 2 , n 3 , -1, n 5 - 1, n 6 + 1) 
+(ni - 1) ni n & f(ni - 2, n 2 , n 3 , -1, n 5 + 1, n 6 - 1) 
-(n 2 - 1) n 2 n 5 f(ni, n 2 - 2, n 3 , -1, n 5 - 1, n 6 + 1) 
-(n 3 - 1) n 3 n 6 /(ni, n 2 , n 3 - 2, -1, n 5 + 1, n 6 - 1) 
+(ni - n 2 - n 3 )n 5 n 6 /(n 1 ,n 2 ,n 3 , -l,n 5 - l,n 6 - 1) 
+n 5 n 6 f(m, n 2 , n 3 , 0, n 5 - 1, n 6 - 1) »i 
-(ni - 1) ni /(ni - 2, n 2 , n 3 , -1, n 5 , n 6 + 1) u> 2 
+ (n 2 - 1) n 2 /(ni, n 2 - 2, n 3 , -1, n 5 , n 6 + 1) u> 2 
-(ni -n 2 - n 3 ) n 6 f(m, n 2 , n 3 , -1, n 5 , n 6 - 1) w 2 
-n% f(ni, n 2 , n 3 , 0, n 5 , n 6 - 1) w x w 2 
-(ni - 1) ni f(ni - 2, n 2 , n 3 , -1, n 5 + 1, n 6 ) u> 3 
+ (n 3 - 1) n 3 /(ni, n 2 , n 3 - 2, -1, n 5 + 1, n 6 ) u> 3 
-(ni - n 2 - n 3 ) n 5 f(n u n 2 , n 3 , -1, n 5 - 1, n 6 ) u> 3 
-«5 /(ni, n 2 , n 3 , 0, n 5 - 1, n 6 ) t«i u> 3 

n 2 , n 3 , 0, n 5 , n 6 ) w 1 w 2 w 3 
+n 6 5(n 3 ) T(n 5 - 1, n 6 - 1, ni + n 2 - 1, w 1 + w 2 , w 3 , 0) 
+n 5 <5(n 2 ) T(n 6 - 1, n 5 - 1, ni + n 3 - 1, iui + w 3 , w 2 , 0) 
-n 5 S(ni) T(n 5 + n 6 - 1, -1, n 2 + n 3 - 1, w 2 + w 3 , Wi, 0) 
-n 6 5(ni) T(n 5 + n 6 - 1, -1, n 2 + n 3 - 1, w 2 + w 3 , wi, 0) 
-8(n 2 ) r(n 6 , n 5 - 1, m + n 3 - 1, iui + u> 3 , u> 2 , 0) u> 2 
+5(ni) T(n 5 + n 6 , -l,n 2 + n 3 - 1, u> 2 + u> 3 , 0) u> 2 
-5(n 3 ) r(n 5 ,n 6 - l,m + n 2 - l,Wi + w 2 ,w 3 ,0)w 3 
+5(ni) T(n 5 + n 6 , -1, n 2 + n 3 - 1, w 2 + iu 3 , wj b 0) w 3 ] . (24) 

This general formula does not work in the case n x = n 2 + n 3 . In the second step we use integration 
by part identities in the coordinate space to fill this hole. We limit ourselves only to a special case 
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of these identities in the form 



= id(i) = J d 3 n J d 3 r 2 J d 3 r 3 (gV\h - h V 2 l9 ) , (25) 



where 



_ _— wi n— W2 r%— wzrz „rM— 1 r ns— 1 r n$— 1 
9 — e '1 '2 '3 J 

^ = rg-vsrvjr 1 . (26) 

The identities id(2) and id(3) 

/(ni,n 2 ,n 3 ,-l,7i5,n 6 ) = [(ni - 1) (ni + n 3 - 1) f(n x - 2, n 2 , n 3 , -1, n 5 , n 6 ) 

-(ni - 1) (n 3 - 1) /(ni - 2, n 2 + 2, n 3 - 2, -1, n 5 , n 6 ) 
+ (n 3 - 1) (ni + n 3 - 1) /(ni, n 2 , n 3 -2,-1, n 5 , n 6 ) 
-(n 5 - 1) n 5 /(m, n 2 , n 3 , -1, n 5 - 2, n 6 ) 
+2 n 5 /(m, n 2 , n 3 , -1, n 5 - 1, n 6 ) w 2 

+5(n 5 )T(n l + n 6 - l,n 3 - 2, n 2 , u> 3 , wi, 0)]/wl , (27) 
/(ni,n 2 ,n 3 , -l,n 5 ,n 6 ) = [-(ni - 1) (n 2 - 1) f(m - 2,n 2 - 2,n 3 + 2, -l,n 5 ,n 6 ) 

+ (n 1 - 1) (m + n 2 - 1) /(ni - 2, n 2 , n 3 , -1, n 5 , n 6 ) 
+ (n 2 - 1) (ni + ri2-l) f{ni, n 2 - 2, n 3 , -1, n 5 , n 6 ) 
-(n 6 - 1) n 6 /(m, n 2 , n 3 , -1, n 5 , n 6 - 2) 
+2 n 6 /(m, n 2 , n 3 , -1, n 5 , n 6 - 1) w 3 

+5(n 6 ) r(n 2 - 2, ni +n 5 - l,n 3 ,wi,w 2 , 0)]/w| , (28) 



replace the main recursion in Eq. (l2*?b for the case n\ = n 2 + n 3 and can be used also for all other 
rii under conditions that n\ > 0, n 3 > or ni > 0, n 2 > 0, respectively. 



V. SUMMARY 

We have demonstrated the advantages of the analytic approach to three-electron Hylleraas inte- 
grals by the calculation of nonrelativistic energy of the ground state lithium atom and the leading 
relativistic corrections. The achieved accuracy is the best to date and this is mainly due to the 
use of much larger basis sets. In fact it is possible to perform calculation with f2 > 12 by using 
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sextuple precision arithmetics. The typical evaluation time in sextuple precision for = 12 is 24 
hours on 2.4 GHz Opteron, and most of the time is devoted to LU decomposition. 

Having precise wave functions, we have calculated leading relativistic corrections and the re- 
sults only partially agree with that of Yan and Drake iflSfl and of King [17]. We are now in position 
to calculate higher order, namely m a 6 relativistic and QED corrections, for example to the lithium 
ground state hyperfine splitting J^lll . However, this involves more complicated Hylleraas integrals 
containing two factors among 1 / rf and 1 / r?- , which have not yet been worked out by the recursion 
methods of the authors. 

Even more interesting is the possible extension of this analytic method to beryllium and 
beryllium-like ions, the 4-electron systems. The use of large Hylleraas basis set will allow for 
a high precision calculation of the wave function, energies and transition rates. For example, 



knowing the isotope shifts, one can obtain charge radii as for the lithium is otop e 12211 . General 
Hylleraas integrals for 4-electron systems has not yet been worked out I23I I2411 . The so called 
double linked basis set, the functions with at most two odd powers of have been used by Biisse 
et al in 12511 to obtain an accurate nonrelativistic energy, but still less accurate than the result of 
Komasa in 11(11 . It has not yet been attempted to calculate relativistic corrections with Hylleraas 
functions as they involve even more difficult integrals. We think, the integration by part technique, 
should allow for the derivation of compact formulas for all 4-electron Hylleraas integrals. 

Our primary motivation for developing Hylleraas basis set is the calculation of higher order 
relativistic and QED effects, and to demonstrate that standard techniques used in relativistic quan- 
tum chemistry, which are based on the multi-electron Dirac-Coulomb Hamiltonian are not correct 
for principal reasons. This Hamiltonian does not not include properly negative energy states. The 
correct treatment has to be based on quantum electrodynamics and several very accurate results 



HQ. 



for few electron ions have already been obtained within the so called \/Z expansion [26, 
Nevertheless, there is no yet formalism which allows for systematic inclusion of negative energy 
states and QED effects for many electron atoms. 
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APPENDIX A: QUADRATURE WITH LOGARITHMIC END-POINT SINGULARITY 

Consider the integral 

1= f dx \Wx{x) + ln(x) W 2 {x)] , (Al) 
Jo 

where Wi are arbitrary polynomials of maximal degree n — 1. We would like to find n nodes Xi 
and n weights Wi such that 

n 

I = J2 W * l W ^) + H*i) W 2 {x i )] ■ (A2) 
i=i 

In general it is a difficult numerical problem to find a solution of corresponding 2 n nonlinear 
equations with j = 1, n 

f 1 1 n 

j dx x^' 1 — — — / Wi x( , (A3) 

r 1 i n 

/ dxx^ 1 lnx = — — = Wix{ lnxj . (A4) 
Jo J i=1 

The work H29I1 solves this problem and proves that Wi are all positive. The solution is as follows. 

One defines 2 n functions <pi 

<pk(x) = x fc_1 , for k = l,n (A5) 
<f>k{x) = £ fc_1 In a;, for k — n + 1, 2n . (A6) 

Consider n points X{ which are not necessarily the solution of Equations (IA3IA4I) but are close to 
them, and construct another set of functions <Xj, t\%, for % = 1, n 



^atijfaix), (A7) 



2 n 

(T ,(■>') 

3=1 
2n 

th{x) =^2p ij <l> j (x) , (A8) 
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such that 



<Ti(Xk) = 0, 

a 'i( x k) = S ik , 
Vi( x k) = Sik, 

VIM = 0. (A9) 
The set of conditions (IA9I) uniquely determines the matrices and If Xk are nodes, then 

dxcri(x) = 0, 

i 

dxrji(x) = Wi . (A10) 

o 



If x k are not exactly the nodes, but are sufficiently close, then according to work [ 29], the iteration 

Xi > Xi 

x i = Xi + J\ t ±-L ) (All) 

J Q dxr]i(x) 

converges to nodes, the solution of Eqs. (IA3IA4I) . The only problem now, is to find a sufficiently 
good initial values for X{. For this one constructs a homotopy t) such that 

(J)k(x,t) = x k ^ 1 for k = 1, n , 

4> k (x, t) = (1 - t) y/x + t x k ' x - n ln(x) for k = n + 1, 2 n . (A12) 

At t = 0, 0) are polynomials in y/x, therefore one obtains = y| where ?/i are nodes for 
Gauss-Legendre quadrature. By slowly changing t from one finds the solution at t = 1. In 
the actual numerical calculations we found that the steps ti = i/100 were sufficiently small for 
the above iteration to converge. This generalized Gaussian quadrature can also be constructed for 
other types of functions including various, even nonintegrable singularities. 
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